%% Specifying the conditions
% Comments and changes made by Selina in May 2020
% Get the current working directory
%based on Timmer, Marcel P., Abdul Azeez Erumban, Bart Los, Robert Stehrer, and Gaaitzen J. de Vries.
%  “Slicing Up Global Value Chains.” Journal of Economic Perspectives 28, no. 2 (May 2014): 99–118. https://doi.org/10.1257/jep.28.2.99
%data from them as well
workDir = pwd;

disp(['Current directory: ', workDir])

% Find the position of the string "jpe" in the workDir
jpePos = strfind(workDir, 'jpe');

% Shorten the workDir to include only up to and including "jpe"
if ~isempty(jpePos)
    workDir = workDir(1:jpePos + length('jpe') - 1);
else
    error('The string "jpe" was not found in the current working directory.');
end

disp(['Shortened directory: ', workDir])

q = 'shares'; 
fdc = 'ALL'; 
var = 'TOT'; 
agg = fullfile(workDir, 'datasets', 'import', 'agg_foreign_shares_modified.xlsx')
pre_def='play'; 
post_def='selina'; 
store = 'y'; % define 'y' if you want to store the final result to an excel file or 'n' if not
OUT = fullfile(workDir, 'datasets', 'import'); % specify the folder where you want to store the results
info = 'foreign shares mfg & non-mfg'; % add any info between '' that you want to give to the exported excel file such that you can recognize it easily
dataDir = fullfile(workDir, 'datasets', 'import');
% load data
for y = 1995:1:2011 % all years are considered (loop runs until the end of this program)

    
c = 41; % number of countries Selina: Last country in WIOT is Rest of the World (Row) - also used
n = 35; % number of industries
fd = 5; % number of final demand categories


%Selina: Those file directly downloaded from WIOT webpage
clear WIOT95_00 WIOT01_05 WIOT06_09

if y==1995 || y ==1996 || y==1997 || y==1998 || y==1999 || y==2000     
    load(fullfile(dataDir, 'WIOT95_00')); % WIOTs  
    %Selina: (row, column)
    A = WIOT95_00(1:n*c,((y-1995)*1641)+1:((y-1995)*1641)+n*c); %intermediate-input block for year y
    D = WIOT95_00(1:n*c,(((y-1995)*1641)+1+(n*c)):((y-1995)*1641)+(n*c)+(5*c)); % final demand block for year y
    x = WIOT95_00(1:n*c,(((y-1995)*1641)+1+(n*c)+(5*c))); % column vector of domestic industry gross output for year y
    va = WIOT95_00(6+n*c,((y-1995)*1641)+1:((y-1995)*1641)+n*c)'; % column (ROW?) vector of domestic industry value added for year y
end

if y==2001 || y ==2002 || y==2003 || y==2004 || y==2005      
    load(fullfile(dataDir, 'WIOT01_05')); % WIOTs  
    A = WIOT01_05(1:n*c,((y-2001)*1641)+1:((y-2001)*1641)+n*c); %intermediate-input block for year y
    D = WIOT01_05(1:n*c,(((y-2001)*1641)+1+(n*c)):((y-2001)*1641)+(n*c)+(5*c)); % final demand block for year y
    x = WIOT01_05(1:n*c,(((y-2001)*1641)+1+(n*c)+(5*c))); % column vector of domestic industry gross output for year y
    va = WIOT01_05(6+n*c,((y-2001)*1641)+1:((y-2001)*1641)+n*c)'; % column vector of domestic industry value added for year y
end

if y==2006 || y==2007 
    load(fullfile(dataDir, 'WIOT06_09')); % WIOTs  
    A = WIOT06_09(1:n*c,((y-2006)*1641)+1:((y-2006)*1641)+n*c); %intermediate-input block for year y
    D = WIOT06_09(1:n*c,(((y-2006)*1641)+1+(n*c)):((y-2006)*1641)+(n*c)+(5*c)); % final demand block for year y
    x = WIOT06_09(1:n*c,(((y-2006)*1641)+1+(n*c)+(5*c))); % column vector of domestic industry gross output for year y
    va = WIOT06_09(6+n*c,((y-2006)*1641)+1:((y-2006)*1641)+n*c)'; % column vector of domestic industry value added for year y
end

if y==2008 || y==2009 || y==2010 || y==2011
    load(fullfile(dataDir, 'WIOT08_11')); % WIOTs  
    A = WIOT08_11(1:n*c,((y-2008)*1641)+1:((y-2008)*1641)+n*c); %intermediate-input block for year y
    D = WIOT08_11(1:n*c,(((y-2008)*1641)+1+(n*c)):((y-2008)*1641)+(n*c)+(5*c)); % final demand block for year y
    x = WIOT08_11(1:n*c,(((y-2008)*1641)+1+(n*c)+(5*c))); % column vector of domestic industry gross output for year y
    va = WIOT08_11(6+n*c,((y-2008)*1641)+1:((y-2008)*1641)+n*c)'; % column vector of domestic industry value added for year y
end

fd_v = D * ones(205,1); %Selina: generates the sum of the final demand across the rows (205 columns) for each country and industry
coef = A*invd(x); % intermediate input coefficients = Selina: how much intermediates are needed to produce a unit of output of a given product ; column of intermediate-input of industry in country over total output (output at basic prices) in industry in country, "Technologiematrix: https://www.youtube.com/watch?v=YdFON0-Je6U - wieviel konsumiert �ber wieviel produziert = wieviel technolgie man braucht"   invd is a package from them - can be found in the source folder
L = inv(eye(1435)-coef); % leontief inverse  %Selina: eye=Identity matrix - rest how one computes the leontief inverse 

% for shares decomposition
if strcmp (q,'shares'); %Selina: compares two strings
    if strcmp (var,'TOT');
        FINAL = diag(div(va,x))* L * diag(fd_v);  %Selina: div is a package from them - divides va by x; is in the online appendix k=F(I-B)^-1 *z
        %Selina: A typical element in k indicates the value added in the
        %production of the first final good by the particular production
        %factor (here all production factors together) employed in country
        %i and industry j
    end   
end



% pre-multiplication
if strcmp (pre_def,'play'); % decide which country and industry groupings along the rows to load - Selina: only this used
    [pre_c] = xlsread(agg,'pre','c3:c43');
    [pre_i] = xlsread(agg,'pre','m3:m37');
end

c_pre_groupings = unique(pre_c(:,1)); % distilling the number of country groupings
i_pre_groupings = unique(pre_i(:,1)); % distilling the number of industry groupings
pre = zeros(size(c_pre_groupings,1)*size(i_pre_groupings,1),1435); % pre-specifying a matrix that has dimension country*industry by 1435

% filling the pre-vector or pre-matrix
for c=1:41 % all countries (includes RoW)
for i=1:35 % all industries
for cg = 1:size(c_pre_groupings,1) % loop across number of country groupings
for ig = 1:size(i_pre_groupings,1) % loop across number of industry groupings
    if pre_c(c)==(c_pre_groupings(cg,1)) && pre_i(i)==(i_pre_groupings(ig,1)) % check whether country and industry grouping match
        pre(((cg-1)*(size(i_pre_groupings,1)))+ig,(c-1)*35+i)=1; % replace pre-specified matrix with a 1 if country and industry grouping match                     
    end
end
end
end
end

% post- multiplication
% Selina: New use this to get mfg and non-mfg industries
if strcmp (post_def,'selina'); % decide which country and industry groupings along the columns to load
[post_c] = xlsread(agg,'post','d3:d43');
[post_i] = xlsread(agg,'post','n3:n37');
end


c_post_groupings = unique(post_c(:,1)); % distilling the number of country groupings
i_post_groupings = unique(post_i(:,1)); % distilling the number of industry groupings
post = zeros(1435,size(c_post_groupings,1)*size(i_post_groupings,1)); % pre-specifying a matrix that has dimension1435 by country*industry

% filling the post-vector or post-matrix
for c=1:41 % all countries (includes RoW)
for i=1:35 % all industries
for cg = 1:size(c_post_groupings,1) % loop across number of country groupings
for ig = 1:size(i_post_groupings,1) % loop across number of industry groupings
    if post_c(c)==(c_post_groupings(cg,1)) && post_i(i)==(i_post_groupings(ig,1)) % check whether country and industry grouping match
    post((c-1)*35+i,((cg-1)*(size(i_post_groupings,1)))+ig)=1; % replace pre-specified matrix with a 1 if country and industry grouping match
    end
end
end
end
end

% Generating the required aggregation

RESULT = (pre * FINAL * post); % pre and post-multiplication of the pre-specified global value chain matrix for aggregation

test = ones(41,41);
for i = 1:41
    test(i,i)=0;
end
test = kron(test,ones(2,1))'; % for getting rid of the domestic part %Selina: Kronecker tensor product between two matrixes, changed from 15 to 2 by Selina


VA_for = RESULT .* test; %Selina: Element-wise multiplication 

sumVA = sum(RESULT,1); %Selina: take the sum 
sumVA_for = sum(VA_for, 1); 

ratio =sumVA_for ./ sumVA;

f_share((y-1994),:)= ratio; %%Selina: same as ratio but but in right column for the year (1995 at top, then 1996, etc.)


end

share = f_share';

if strcmp (store,'y');  % storing country results, selected industry
    filename = fullfile(OUT, [info, '.xlsx']); 
    [num format] = xlsread(agg,'format','a1:s84');
    clear num;
    xlswrite(filename,format,'share','a1');
    xlswrite(filename,share,'share','c2'); % save result in sheetname that can be traced back to the criteria  

end

exit(0);

